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ABSTRACT 

The new results of the evolution of the synchrotron peak for Mrk 421 are 
mostly likely related to the particle acceleration process. In order to account 
for the above results, we present a model of blazar variability during the flare 
in which the emission comes from accelerating electrons. A diffusion advection 
equation of the electron energy distribution is derived to calculate the spectrum 
and light curve of synchrotron radiation. We present that the observed shifts 
of the synchrotron peak moving to higher energies during the flare are caused 
by shock acceleration. The observed relation between changes in the fluxes at 
specifled frequency ranges and shifts of the peak position is fltted to constrain 
the physical parameters of the dissipation region. 

Subject headings: accleration of particles-galaxies: active- BL Lacertae objects: 
indivual(Mkn 421)-Plasma:turbulence-radiation mechanisms: nonthermal- 
shock waves-X-rays: general 
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1. Introduction 

Among active galactic nuclei, blazars are characterized by a high and variable degree 
of polarization and a flux variability often occurring on very short time scales. The 
observed continuum is dominated by nonthermal emission, where the emission deriving 
from synchrotron and self-Compton is enhanced by relativistic beaming (Blandford & Rees 
1978). Some blazars exhibit not only intensity variation extending from radio waves to 
gamma rays, but also spectral variation as a function of the flux level in many frequency 
bands. 

An interpretation of the spectral variations can be envisaged in an inhomogeneous 
jet scenario of the kind proposed by Ghisellini et al. (1985) and Celotti et al. (1991). In 
this model the structure of the jet is considered, the magnetic fleld and the relativistic 
electron density are assumed to decrease with distance, and the maximum synchrotron 
frequency is assumed to decrease with distance, accounting for a variability time shorter at 
larger synchrotron frequencies. The predicted variability pattern is caused by a disturbance 
(e.g. shock) traveling down the jet. This perturbation is assumed to produce only a flxed 
enhancement in magnetic field and particle density and to unchange the shape of the 
electron distribution. The overall spectrum is the superposition of the located spectra 
emitted by each shce of the jet. As it has been apphed to BL Lac objects, the model does 
not consider the evolution of the relativistic electron distribution affected by radiation 
losses, new injection of relativistic electrons and shock waves. However, simultaneous 
observations of blazars over a wide spectral range suggest that significant radiation losses 
and injection of relativistic electrons occur in fiares. For example, the variability of 
PKS2 155-304 at X-ray, EUV and UV bands is decreased in amplitude and with significant 
delays approximately satisfying a, relation (Urry et al. 1997). In particular, the 
radiation loss time for electrons that emit optical and X-ray emission in blazars is probably 
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far less than the hght (or the shock) crossing time of the emitting region. Consequently, 
if the electron distributions responsible for high energy emission are kept by the shock 
acceleration processes, then radiation losses and the injection of electrons occur in the 
vicinity of the shock during the flare and must be considered in the theory. Many models 
to reproduce the spectral variabihty have been developed ( Mastichiadis & Kirk 1997; 
Georganopoulous & Marscher 1998; Kirk, Rieger & Mastichiadis 1998; Wang ct al. 1999; 
Ghiaberge & Ghisellini 1999; Li & Kusunose 2000). The observed spectrum and time delays 
between the light curves at fluxes at different frequencies are believed to be produced by the 
electron distribution at different stages of evolution after episodic electron injection phases. 
However these models can not account for a synchrotron peak drifting to higher energies 
at X-ray wavelengths during the rising phase of the flare in Mrk 421 (Fossati et al. 2000). 
The above hard lag is most likely related to the particle acceleration process. The role of 
particle acceleration by shock waves has been considered by Kirk ct al. (1998). In this 
model, electron acceleration and radiation zones are separated. Electrons are continuously 
accelerated at the shock front, and subsequently drift away from it in the downstream and 
emit most of the radiation. The emission from the acceleration zone is ignored. 

Through shocks can quickly accelerate particles to very high energies, this requires 
the existence of some scattering agent to force repeated passage of the particles across 
the shock. The most likely agent for scattering is plasma turbulence or plasma waves. 
However, the plasma turbulence needed for the scattering can not only accelerate particles 
stochastically (second order Fermi acceleration), but can cause particle diffusion in the 
downstream zone quickly. If the dissipation of the bulk energy of blazar jet into particles 
includes shocks and plasma turbulence, it seems difficult to distinguish acceleration zone 
and emission region. The emitting particles in downstream could be re-accelerated by 
shock waves under the scattering of plasma turbulence. Therefore we try to consider the 
synchrotron emission from a single dissipation region which includes shock and turbulence 
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processes and radiation losses. 

In Section II we present a stochastic differential equation to describe a particle 
acceleration and cooling processes, together with the assumptions made. Then we derive a 
diffusion-advection equation of particle energy distribution which simulates the temporal 
evolution of the particle radiation occurring when shock waves and plasma turbulence are 
produced in a dissipation region of a jet. In section III we apply our model to Mrk 421 to 
explain the hard lag observed in the X-ray light curves. Finally, we draw our conclusions in 
section IV. 

2. The Model 

We focus on the particle acceleration and cooling processes and assume that pitch-angle 
scattering maintains an almost isotropic particle distribution. We also assume that the 
shock in a relativistic blazar jet is nonrelativistic, and assume that the light crossing time 
of the dissipation region is short compared with the intrinsic cooling and acceleration 
timescales. It is shown that the observed variability is determined by the intrinsic cooling 
and acceleration processes. We use a homogeneous model in which both the magnetic field 
and the particle distribution function are assumed homogeneous through the dissipation 
region, and consider only the synchrotron radiation of the accelerating particles, leaving the 
more involved computation of inverse Compton emission to future work. 

2.1. Evolution of particle energy distribution 

Firstly we use a stochastic differential equation (SDE) to describe a particle acceleration 
and cooling processes, as has been used by some authors (Kriills & Achterberg 1994; 
Mastichiadis & Kirk 1997). Then we derive a diffusion-advection equation of particle energy 
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distribution (or Fokker-Plank equation). Because the particle distribution is assumed to 
be isotropic in space and momentum, we consider only the evolution of the particle energy 
given by 



de 
'dt 



''-gam OSS 



The energy losses Sioss is mainly determined by synchrotron radiation and is a deterministic 
process in homogeneous model. It has the form 



The energy gains e^am is due to shock wave and plasma turbulence, and is assumed to have 
the form 



where a is a stochastic variable which is determined by mean value (a) = a and a variable 
given by the Gaussian noise process, namely it defined by the value {({t)) = and the 
autocorrelation function (C(^)C(^')) — '2D5{t — t'), where the coefficient D denotes the noise 
intensity which is assumed to represent the stochastic influence of plasma turbulence. This 
means that the turbulence will influence the rate of shock acceleration stochastically and 
cause the particle energy gain to be a random process. To that end one writes Eq(2) as an 
SDE 



_^«£-/3£^ + £CW- (4) 

For simplify we have assumed that a, {3 and D change slowly compared to particle cooling 
and acceleration processes, and treat them as constants. 
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The Fokker-Plank equation corresponding to the system of SDE given in Eq(4) are 
(Stratonovich 1967) 



df{e,t) 



dt 



d_ 
de 



+ D 



92 



(5) 



{ae - (3e^ + De)f{e,t) 

where f{£,t) denotes a conditional probabihty density p{e,t/eo,to) for any initial energy 
£0 at ^0- H the number A^o of particles undergo acceleration, the evolution equation of the 
particle number density N{e,t) = NQf{e,t) is similar to Eq(5). In order to include particle 
escape and injection, we are easy to extend the evolution equation as (Schlickeiser 1984) 



dN(e,t) d r, 2x,.. xi r^d'^ [2,./ ^^ N(e,t) , 

-^f^ = [(«^ - f^^')N{e, t)] + D— [e'N{e, t)\ - -jM + Q{e, t). (6) 

We consider only the intrinsic cooling and acceleration processes and assume that the 
timescale of particle is longer than that of the intrinsic physical processes. We ignore the 
particle escape term in Eq(6). Otherwise we assume that the particles in the dissipation 
region are accelerated from the initial time to and no new particle is injected in the following 
time. The injection term in Eq(6) is also ignored. We scale the energy s and the time t 
as e = EiX and t = a~^T, where the critical energy ei = aP~^, and also define k, = Da^^ 
which is the ratio of shock acceleration timescale to turbulent acceleration timescale. We 
can now rewrite the evolution equation of the particle distribution as 



dN{X, t) 
97 



d 



dX 



rj2 

X - + 5X] N{X, t) + [X^N{X, t) 



(7) 



We now turn to finding the solution of Eq(7) describing the evolution of 
particle distribution. Assume the initial energy distribution of the particles to be 
N{X,0) = S{X — xq). We firstly consider the solution of Eq(7) in the case of weak 
interaction of plasma turbulence k <^ 1. 
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We write X &s X ^ x{t) + where x{t) = {X) = J XN{X,T)dX. Multiplying 
Eq(7) by X and X'^ and integrating X, we obtain the equation of moments 



dx{T) 
dr 



x{t) -X {t), 



da{T) 
dr 



2[l-x(r)]a(r) + 2a;2(r), 



(9) 



where (j(r) = (y^) = / y'^N{X,T)dX. Other high order moments can be ignored due to 
K <^ 1. Therefore, the solution of Eq(7) is approximately given by a Gaussian distribution 



N{X,t) 



: exp 



[X ,ir)Y 



(10) 



where x{t) and cr(r) are respectively the solutions of Eq(8) and Eq(9) which are given by 



l)e- 



and 



a{T)^C[x{T)]-C{xo). (12) 

The function C{x) is defined as 

C{x) = 2{x - x^f In + + " ^) (13) 

shown in Fig.l. Clearly the energy distribution is narrow, centered at a typical energy 
Six{t) and with width Si^ nair), quickly evolving with time to £/ and Swfn^ where Si 
corresponds to the highest energies. It is shown that the particle acceleration is limited by 
radiation losses 
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We now consider the solution of Eq(7) for the general value of k. It is noted that Eq(7) 
has a stationary solution given by 

iV,(X)=CX^-iexp(-^), (14) 

where C is a normalization constant. The function Nc{X) has a peak at X — 1 — k and a 
narrow shape with decreasing k, which is shown in Figure 2. It indicates that the turbulence 
leads the spread of the particle energy distribution. The formation of a peak is caused by 
the interplay of shock wave acceleration, turbulence acceleration and radiation losses 

The time-dependent solution can be constructed by eigenfunctions and eigenvalues of 
Eq(7) (Schenzle & Brand 1979): 

N{X,t)^Y.^-'^^C,^P,^{X), (15) 

where Px^^lX) denotes the eigenfuctions of the operator A{X) with eigenvalues A^, Cx^ are 
expansion coefficients, 



A{X)Px(X) 



d 



dX 



{X - X^ + 5X)PxAX) 



+ K 



dx^ 



X'Px,{X) =X,Px,iX), (16) 



where A{X) is not an Hcrmtean operator. The condition of square-integrability and the 
correlation for the eigenfunctions are 



N,{X)[Px^{X)] dX^l, Py^{X)Px^{X)dX^S{X'^-X,), (17) 

where the function Px'^{X) = Px^{X) / Ns{X) . From the above condition, the allowed 
eigenvalues and eigenfunctions are given by 
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Xn — n — Kri^, (18) 

Px„{X)^Ns{X)X-^M{-n,-2n + l + K-\XK-^) (19) 

for the discrete part of the eigenvalue spectrum, where n is the positive integer subject to 
the condition n < and a continuous spectrum for > ^ with the corresponding 

eigenfunctions 



Px,iX) = NsiX)X"'-^^Uiia - — , 2ia + 1, Xk'^), (20) 

where M{a, b, X) and U {a, b, X) are Rummer's function (Abramwitz & Stegun 1970) and 
a = [Xi,5k - ^K-') . The expansion coefficients Cx are determined from the initial 
condition 



roo 

^'^^Jo N7\X)Px,imiX,0)dX. (21) 

One finally finds 

/oo 
, e-'^^C,^P,^{X)dX,. (22) 

We now turn to consider the evolution of the particle distribution when the particle 
acceleration stops and ignore the particle escape. The evolution equation of the particle 
distribution is given by 



dN{e,t) _ d 
di ~ ds 



Pe^N{e,t) +Ns{e)d{t-t,), 



(23) 



where Ns{£) and ts are respectively the particle distribution and the time when the 
acceleration stops. This case corresponds to the decay phase. The solution of Eq(23) is 
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N{e,t) = Ns{eb-^)b 



b = [l- I3e{t - ts)] > 



(24) 



which is an extension of the solution given by Kazanas et al. (1998). 



3. Application to Mrk 421 



A peaked spectrum from observations constrains particle energy distribution to be 
narrow. It indicates that shock acceleration will dominate particle acceleration processes 
corresponding to -C 1. In the following text we focus on the calculation of light curves 
and spectra in the case of k <^ 1. We assume that the initial energy distribution of electrons 
to be single energy distribution with initial energy xqEi. The synchrotron emission is 
calculated from the energy distribution of electrons which is parameterized by the critical 
energy Si. We first estimate the parameter ei from the timescales of shock acceleration and 
radiation losses. The timescale of shock acceleration can be estimated in terms of the 
light crossing timescale At. We assume the timescale a^^ to be ( times the light crossing 
timescale, e.g., = C^i, one obtain 



We now calculate the synchrotron emissivity as a function of time and frequency 
F^{u,t) with 




(25) 



The critical frequency is given by 




(26) 
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F^{u,t) oc J deN{e,t)P{u,e). (27) 

P(i/, e) is the single particle synchrotron emissivity which has a approximate function given 
by Kaplan & Tsytovich (1973): 

P{u,e)^a,V3l^iy Fl^£j, (28) 

where = 27rV3e^c~^ULsm9 is a constant and Uc — Se^ULsm9/{2mlc^), with 

vl — eB/{27rmec) the electron Lamor frequency and 9 the angle between the magnetic field 

direction and the line of sight. The function F(X) is F(X) = X^/^exp(— X). 

Because the particle distribution is a narrow Guassian distribution, it can be 
approximated as a Delta function, e.g., N[e,t) = iVo5(e — Ec), where Sc is the central 
energy of Guassian distribution which is evolving with time {Sc — x{t)si). The synchrotron 
emission is approximated as Fj, oc P^u^Sc)- In Figure 3 we plot the synchrotron spectrum 
vFi, of accelerating electrons at different times. Clearly the synchrotron peak shifts to 
higher frequencies during the rise. It indicates that the presence of a hard lag is caused by 
the emission from the electrons which are being accelerated to higher energies by shock 
wave. The highest energy is determined by the acceleration and cooling rates. 

When the shock acceleration stops, the particle distribution will evolve according to 
Eq(24) due to radiation losses. The corresponding synchrotron emission is given by 

F^oc r NJeb-^)b-^P(u,e)de, (29) 
Jo 

where e* is a cut-off frequency which comes from the limit of 6 = [1 — (3e{t — td)] > 0, e.g., 
e* — P~^{t—ts)'^. Using a Delta function to approximate Nd, e.g., N(i{eb~^) — NoS{eb~^—ec), 
we obtain the synchrotron emission as 
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oc r 5{eb-^ - ec)b-^P{u, e)de = P{u, St), (30) 
Jo 

where = e^l + ScP{t — ts)\~^ and Sc is the central energy of particle Guassian distribution 
when the shock acceleration stops. Clearly the peak of vF^, decreases to lower frequencies 
as soon as the shock is over. 

The brightness of synchrotron emission at specified frequency ranges (1^1,1^2) for the 
flare is simply given by 

Lf = J F^dp (X J P{p,ec)dp. (31) 

The peak frequency of synchrotron spectrum i/F^ locates at i/p — ^v^.. The relation of Vp 
and ^/ is 

Lf oc Q{up, U1U2) = Up X^e'^dX, (32) 

where the function Q{i'p, ^1^2) is shown in Figure 4. The peak frequencies show an 
obviously linear relation with the fluxes in 0.1-lOkeV ranges. This relation reveals the 
formation of a narrow particle energy distribution by shock wave and the subsequent 
evolution by radiation losses. 

In order to show the above results for application, we thus study the relation between 
changes in the brightness and shifts of the peak position during the flare based on data of 
Mrk 421 on 1998 April 21 (Fossati et al. 2000). We introduce a steady spectrum 6(z/) which 
does not take part in the flare. The observed brightness Lq^s will include two contributions 
of the flaring and steady components given by 



Lobs = A [Q{^p, ^1,^2) + K{vi, i/2)] , 



(33) 
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where K{i/i, 1/2) — j^l h{i')di>, and A — agV^ (|j Uq (^j c/St which is given by 

^ - X 1°-'° il^) (oij) (34) 
where Uq is the number density of accelerated particles. B and it! are the magnetic field 
and size of dissipation region respectively. Dj^ is the luminosity distant of source to the 
observer. We take the Hubble constant to be 50A;ms^^Mpc^^and the redshift of Mrk 421 to 
be 0.031. The fits of the model are shown in Figure 5. The values of the model parameters 
are A = 10.0 x IQ-^^ergs-'^cm-'^keV-^Qjid K = 0.08, 0.23, OAlkeV for the 0.2-1, 2-10, 
0.1-lOkeV fluxes. Thus the relation between the spectral evolution and the flux variability 
during the flare can constrain the physical parameters of dissipation regions. It is important 
to notice that the inclusion of a steady component is crucial to flt the above relations. The 
flts of the model indeed achieve the deconvolution of the spectral energy distribution into 
different contributions. The changes of 1^2) with different frequency ranges show that 

the steady emission concentrates at 2-lOkeV ranges. 

The evolution of synchrotron emission with times is given by Eq(30), where the central 
energy St evolves during the acceleration and post-acceleration phase as 

^t = \ , (35) 

[ St = £iy{Ts)[l + y{Ts){T - Ts)]~^ t > ts 

where y(rs) = [1 + (xq^ — l)e~'^]'~^, r — at — ^ and Tg — Figure 6 shows the function 
F{u/uc) of light curves in different frequencies. We flnd the following interesting results. 
During the acceleration, the higher energy emissions lag the lower energy ones. The light 
curve is approximately symmetric in the lower energy bands, and it becomes increasingly 
asymmetric at higher energies. During the post-acceleration, the light curves are traced in 
the opposite way. The higher energy emissions first rise and then decays rapidly. The light 
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curve is symmetric in various energies due to the single rise and decay timescale determined 
by cooling timescale. It should be also noted that the emissions in frequency ranges of 
^i'iXq < u < have the light curves of the peak shape due to the maximum value of 
function F{v/v^ occurring v — + {x^^ — 1)6""^]. The hght curves m. v > vi bands 
successively increase to the maximum value of F{yjv^ until the shock wave acceleration 
stops and then rapidly decay due to particle radiation losses. If the particle escape is not 
ignored, the variability amplitude in lower energies will decreases in the post-acceleration 
(Wang et al. 1999). If there is new particle injection during the shock acceleration, higher 
energy particles will be more than lower energy ones. The larger amplitude variability 
appears at higher energies. 

4. Conclusions 

The recent BeppoSAX observations of Mrk 421 (Fossati et al. 2000) provide important 
information to understand particle acceleration processes. Within a single emission region 
scenario for blazar jets, we have studied the time dependent behavior of the particle 
distribution affected by the particle acceleration and radiation losses, and calculated the 
form of hght curves and spectra at different times. We have presented that the observed 
shifts of the synchrotron peak moving to higher energy during the flare are caused by 
shock acceleration. The accelerating particles follow a narrow Guassian distribution with a 
central energy quickly evolving with time to the highest energy when the shock acceleration 
dominates the turbulent acceleration. The highest energy is limited by radiation losses. Our 
results are important for the observed fast dissipation region of blazar jets where the light 
crossing time of the region is shorter than the particle acceleration and cooling time scales. 
The observed fast variabihty indicates the particle acceleration and coohng processes. 

The observed relation between changes in the fluxes at specifled frequency ranges and 
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shifts of the peak position during the flare can estimate the magnetic field and the number 
density of accelerated particles in the dissipation region. 

An energy dependence of the shape of the light curve observed in Mrk 421 during the 
flare (Fossati et al. 2000) can be connected to the particle acceleration and cooling time 
scales. During the acceleration, the observed light curves are expected to be symmetric at 
low energies where the acceleration time is similar to the cooling time and much longer than 
the light crossing time. The asymmetric light curves (faster rise) occur at higher energies 
where the acceleration time is shorter than the cooling time and is comparable to the light 
crossing time. During the post-acceleration, the light curve is symmetric where the rise and 
decay timescales are equivalent to cooling time scale. 

With the steady state solution of the diffusion- ad vection equation of particle energy 
distribution, we have demonstrated that the turbulent acceleration leads the particle energy 
diffusion. The strength of turbulent acceleration determines the width of energy diffusion 
which modifies the particle distribution with respect to the pure shock acceleration case. 

We acknowledge the support of National Astronomical Observatories grants 
(99-5102CA), Chinese Academy of Sciences (CAS). 
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Fig. 1. — Evolution of the squared variance of particle Guassian energy distribution 
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Fig. 2. — Steady state particle distributions for different turbulent parameters 
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Fig. 3. — Synchrotron spectral evolution as the central energy of particle distribution is 
evolving to higher energies 
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Fig. 4. — Synchrotron peak energy Up versus the fluxes Q{i'p,i'i,i'2) in specifled frequency 
ranges [vi, 1^2] 
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Fig. 5. — Synchrotron peak energy versus the fluxes in specifled frequency ranges for 1998 
datasets of Mrk 421. The sohd hues represent the best fits of the model 
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Fig. 6. — Multiwavelength light curves calculated from F{y/vc) in the case of Xq — 0.1 and 

Ts = 6.0 



